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Many plasmas of interest to the astrophysical and fusion communities are weakly collisional. In 
such plasmas, small scales can develop in the distribution of particle velocities, potentially affecting 
observable quantities such as turbulent fluxes. Consequently, it is necessary to monitor velocity space 
resolution in gyrokinetic simulations. In this paper, we present a set of computationally efficient 
diagnostics for measuring velocity space resolution in gyrokinetic simulations and apply them to 
a range of plasma physics phenomena using the continuum gyrokinetic code GS2. For the cases 
considered here, it is found that the use of a collisionality at or below experimental values allows for 
the resolution of plasma dynamics with relatively few velocity space grid points. Additionally, we 
describe implementation of an adaptive collision frequency which can be used to improve velocity 
space resolution in the collisionless regime, where results are expected to be independent of collision 
frequency. 

I. INTRODUCTION 

Velocity space dynamics are often important in the weakly collisional plasmas prevalent in astrophysics and fusion 
applications, leading to the necessity of a kinetic treatment. Since the kinetic description requires a six-dimensional 
phase space, simulating weakly collisional plasma processes can be computationally challenging. Employing the 
gyrokinetic ordering^ 1 * 2 * 3 ! reduces the dimensionality by eliminating gyrophase dependence, but we are still left with 
a high-dimensional system. Consequently, one would like to know how many grid points are necessary along each 
dimension, particularly in velocity space, in order to resolve a given simulation. 

In the absence of collisions or some other f orm of dissipation, the distribution of particles in velocity space can 
develop arbitrarily small-scale structure d 4 * 5 * 6 * 7 !. This presents a problem for gyrokinetic simulations, as an arbitrarily 
large number of grid points would be necessary to resolve such a system. Of course, all physical systems possess a 
finite collisionality, which sets a lower bound on the size of velocity space structures and, therefore, an upper bound 
on the number of grid points required for resolution. We would like to know how sensitive the plasma dynamics are 
to the magnitude and form of the velocity space dissipation. In particular, we would like answers to the following set 
of questions: Given a fixed number of grid points, how much dissipation is necessary to ensure a resolved simulation? 
Alternatively, given a fixed amount of dissipation, how many grid points are necessary to ensure a resolved simulation? 
Furthermore, what measurable effect, if any, does the addition of dissipation have on collisionless plasma dynamics? 

These questions have been addressed for very few plasma processes^!, in large part due to the computational expense 
involved with such a study. In this paper, we propose computationally efficient diagnostics for monitoring velocity 
space resolution, and we apply these diagnostics to a range of weakly-collisional plasma processes using the continuum 
gyrokinetic code GSf^D. With the aid of these diagnostics, we have implemented an adaptive collision frequency that 
allows us to resolve velocity space dynamics with the approximate minimal necessary physical dissipation. We find 
that the velocity space dynamics for growing modes are well resolved with few velocity space grid points, even in 
the collisionless limit. Including a small amount of collisions (y <C ui) is necessary and often sufficient to adequately 
resolve nonlinear dynamics and the long-time behavior of linearly damped modes. 

The paper is organized as follows. In Sec. II we discuss velocity space dynamics in gyrokinetics and provide examples 
illustrating the development of small-scale structure in collisionless plasmas. Sec. Ill contains a brief overview of 
the GS2 velocity space grid and its dissipation mechanisms. We describe diagnostics for monitoring velocity space 
resolution in Sec. IV and apply them to a number of plasma processes. In Sec. V, we introduce an adaptive collision 
frequency and present numerical results. We discuss our findings in Sec. VI. 
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II. GYROKINETIC VELOCITY SPACE DYNAMICS 



GS2 solves the coupled system consisting of the low-frequency Maxwell's equations and the nonlinear, electromag- 
netic gyrokinetic equation with a model Fokker-Planck collision operator: 
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is the sum of the curvature and V£? drift velocities, 



v^fbxV ( X ) R (4) 

is the generalized E x B velocity (including both the E x B drift and the drift due to the motion of the perturbed 
magnetic field), 

X = $ - V • A (5) 

c 

is the generalized electromagnetic potential, (-) R denotes a gyro-average at fixed guiding center position R, and 

= F M (l ?pj (6) 

is the lowest order expansion of a Maxwell-Boltzmann distribution. The exact form of the collision operator, C[h], 
used in GS2 is discussed briefly in Sec. Ill and described in detail in Refs. ITT1 and [T2l 

We can group the various terms in the gyrokinetic equation (JlJ into three distinct categories: source terms, labeled 
by iS, which typically drive large-scale structures in velocity space; convection terms, labeled by /C, which lead to 
phase-mixing and the development of small-scale structures in velocity space; and dissipation, given by the collision 
operator, which smoothes the distribution function towards a shifted Maxwellian velocity distribution. In general, 
the structure that develops from the balancing of these terms can be quite complicated. However, we can gain insight 
into how small-scale velocity structures develop by considering simplified collisionless systems. 

In the absence of collisions, arbitrarily small scales can develop in velocity space. This is a result of phase-mixing, 
arising due to convection in real space^. As a simple example of this phenomenon, we include in Appendix A a 
calculation of the perturbed distribution function for the collisionless ion acoustic wave in a slab. The result, quoted 
here, illustrates the tendency of collisionless plasma processes to drive small-scale velocity space structures: 

h(z, v {] , t) = e <fe li(*-V)£( W|| ) + n(z, v h t), (7) 

where the overbar on fi indicates an average over perpendicular velocities. The quantities Q and Ti. are explicitly 
derived in Appendix A. Here, it is sufficient to note that both Q and Ti are smooth functions of the parallel velocity. 
The presence of the oscillatory factor e II ii * in the first term (often called the ballistic term) leads to the development 
of a characteristic wavelength in velocity space that decreases inversely with time. The amplitude of this ballistic 
term remains comparable to the second term in Eqn. ([7]) for all time, leading to the development of large amplitude 
oscillations of the distribution function at arbitrarily small-scales in velocity space. A snapshot of this behavior at 



t = 10 (k\\v ty i) 1 is shown in Fig. 1 
The same calculation carried out fi 



for the collisionless ITG mode in a slab yields a distribution function with a similar 
ballistic term component. However, since this mode is linearly unstable, there is also a term describing large-scale 
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FIG. 1: Plot showing development of fine structure in /(ti||) at t = 10 (k\\Vt,i) 1 . The parallel velocity on the horizontal axis 
is normalized by v t h, and f(vu) was initially a Maxwellian. 

structure in velocity space whose amplitude grows in time to dominate the distribution function. As a result, no 
significant small-scale structure develops. This is a typical feature of linearly growing modes in the collisionless limit. 

Of course, all physical systems have a finite collisionality. The dissipation arising from this collisionality is critically 
important: It is a necessary requirement for the existence of a statistically steady state^H and it sets a lower bound 
on the scale-size of structures in velocity spaced A simple estimate for the scale-size of velocity space structures can 
be obtained by assuming a steady state and balancing the collisional term with the other terms in the gyrokinctic 
equation. Noting that C ~ vv\ h d1 (see e.g. Refs. fT3l or IT4|) . we find 

- ~ (8) 
Vth V w 

where v is the collision frequency, u> is the dynamic frequency of interest, v t h = \/2T/m is the thermal velocity, and 
6v is the scale-size of fluctuations in velocity space. This estimate predicts that velocity space structures much smaller 
than the thermal velocity develop in the weakly collisional limit, v < as we would expect from our consideration 
of simplified collisionless systems. 

III. GS2 VELOCITY SPACE 

In order to fully understand the velocity space resolution diagnostics described in later sections, it is necessary for 
the reader to have a basic knowledge of the way in which velocity space dynamics are treated in GS2. To that purpose, 
we now give a brief explanation of the velocity space coordinates and dissipation mechanisms employed in GS2. 
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A. Velocity space coordinates 



Only two velocity space coordinates arc necessary in gyrokinctics because gyroaveraging has eliminated any gy- 
rophase dependence. Fundamentally, GS2 uses energy, E, and a quantity related to magnetic moment, A = fi/E, as 
its velocity space coordinates. This choice eliminates all velocity space derivatives from the collisionless gyrokinetic 
equation and simplifies the discretization of derivatives in the model collision operator. Consequently, the spacing 
of the velocity space grid points is chosen to provide accurate velocity space integrals while satisfying the necessary 
boundary condition at particle bounce points. 



1. Energy grid 



The volume element in velocity space can be written 

B 



dX 



VI - AS, 



dv v 



(9) 



o Jo 



where i? is the gyroangle and a denotes the sign of Vu . Until recently, the energy grid in GS2 followed the treatment 
of Ref. 1151 which places energy integrals in a convenient form by a change of variables to 



X(x) 



zxe 
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where x = v/v t h- This transforms the range of integration from x € [Q, oo) tol£ [0, 1): 

dX rl 



d A v = 
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dX e x 



(11) 



The integration domain is split into the subintervals [0,Xo) and [Xo,l), with the perturbed distribution function 
assumed to be approximately Maxwellian on [Xq, 1). Gauss-Legendre quadrature rules^'are then used to determine 
the location of the grid points in the interval [0, Xq). 

This energy grid provides spectrally accurate energy integrals (i.e. error ~ (1/N) , where N is the number of energy 
grid points), provided the integrand is analytic in X over the integration domain (see e.g. Ref [T7)l . Unfortunately, this 
is seldom the case. To understand why, we consider the functional form of x{X). Taylor expanding X about x = 0, 
we find X ~ x 3 , or equivalently, x ~ X 1 / 3 . This indicates a branch cut in x originating from X = 0, so that most 
functions of x are non-analytic at X = 0. Furthermore, one can show that x — > oo like x ~ -y/ln (1/ (1 — X)) as X — > 1, 
making x non-analytic at both ends of the domain in X . This can be seen in Fig. [2] where we examine x(X). The fact 
that x possesses singularities at the endpoints of the domain in X means that the integration scheme is not spectrally 
accurate for most integrands of interest (especially since the Bessel functions Jo(k±v±/Q) and Ji(k±v±/Q), which 
are non-analytic at X = and X = 1, appear in all integrals of the distribution function at fixed particle position r). 
This is demonstrated in Fig. |3j where we examine the accuracy of the numerical integral of /i(R) = Fm (at fixed r) 
as we vary the number of velocity space grid points. 

In order to achieve spectral accuracy, we have implemented a new energy grid. We begin by splitting the velocity 
integration into two separate integrals: 



dx x 2 G{x) 



dx x 2 G(x) 



dx x 2 G(x), 



(12) 



where xq is a free parameter and G{x) is the function we wish to integrate. On the first interval, (0,xo), we use 
Gauss-Legendre quadrature rules in x to obtain grid locations. Note that use of x as our integration variable ensures 
that the integrand x 2 G(x) will be analytic as long as G is analytic in x over the interval. 



For the interval (a;g,oo) we make the change of variable y 



Xq to transform the integral to 



dx x 2 G(x) = 



1 



dy e y 



z*yJy + xlG(x) 



(13) 



We then use Gauss-Laguerre quadrature rules in y to obt ain grid locations. Note that the volume clement is analytic 



within the domain of integration, as is x{y) 
analytic function of x. 



\/x 2 + y, so that the integrand will be analytic as long as G is an 
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Our use of spectral integration techniques (i.e. Gaussian quadrature), coupled with the analyticity of our integrand 
for well-behaved functions G(x), ensures the spectral accuracy of our integration scheme. While an exponential order 
of convergence is assured, the rate of convergence depends on the exact nature of the integrand and our choice of the 
parameter xq. In general we choose Xq > 2.5 so that the branch cut at y = —Xq is sufficiently far from the domain of 
integration in y to minimally impact the rate of convergence. We demonstrate the spectral accuracy of the scheme 
and determine the rate of convergence in Fig. [3] It is worthwhile to note that for few grid points (< 8 in Fig. [3]) the 
grid given in Ref. [15|may be more accurate. This is because the energy variable X eliminates velocity-dependence of 
the volume clement (when solving for the normalized distribution function h = h/F ), while the new v-space integrals 

described here have the velocity-dependent volume element x e~ x that must be integrated regardless of the form of 
h. 



2. Lambda grid 

For systems with curved magnetic field lines, special care is also required when dealing with >P^. There are two 
reasons for this: the grid points provided by Gaussian quadrature rules are concentrated near the endpoints of the 
domain, whereas one would like them to be concentrated at the trapped-passing boundary; and one must ensure that 
the proper boundary condition (i.e. /(u|| = + ) = /(i>|| = - )) is satisfied at each of the bounce points. Consequently, 
the A-grid is divided into two regions corresponding to trapped and untrapped particles, respectively. 

For values of A such that < A < 1/B max , the corresponding particles are untrapped by the magnetic potential 
well. In this region of velocity space, the integration variable £ = \/l — A £>„,„ ,. is chosen. It is similar to pitch-angle, 
but it has no spatial dependence. Similarly to the energy, Gauss-Legendre quadrature rules are used to obtain the 
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FIG. 3: Plot showing relative error in numerical integral of Jo(k±v±/Qo) as the number of energy grid points is varied. 
The error due to the integration scheme of Ref. [TF| (solid line) obeys a power law in the number of grid points, with an 
exponent of approximately —2.4. The new integration scheme detailed here (dashed line) has error approximately proportional 
to (0.42iV)~ a42JV , where N is the number of energy grid points. Note that the minimum error in our integration scheme 
approaches 1CP 16 , which is a limitation imposed by double precision evaluation of the Bessel function. 

location of grid points in £. This naturally provides a concentration of gridpoints near the trapped-passing boundary. 

For values of A such that 1/B max < A < l/B m i nj the corresponding particles are trapped by the magnetic potential 
well. In the trapped region, grid points are chosen to fall on bounce points in order to allow for the enforcement of 
boundary conditions. Mathematically, this means that for each value of 9, there must be a corresponding A such that 

V J9) , 

£(0) = J±l = y/l- XB (9) = 0, (14) 

V 

where 9 gives the position along the unperturbed magnetic field line and £ is the pitch-angle. This choice of A values 
also leads to a concentration of grid points near the trapped-passing boundary. A typical GS2 grid layout for a system 
with trapped particles is shown in Fig. [4] It should be noted that the A integrals, like the energy integrals, are 
spectrally accurate, provided the distribution function is analytic in £. 

B. Velocity space dissipation 

Some form of dissipation is often necessary to prevent the formation of arbitrarily small-scale structures in velocity 
space. This can be achieved either through artificial numerical dissipation or through implementation of a model 
collision operator. Both options are available in GS2. 
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FIG. 4: Typical velocity space grid used in GS2. Grid points are concentrated near the trapped-passing boundary (whose 
location varies with 9) and at lower energy values where the Maxwellian weighting dominates. Red (blue) grid points are 
sample A (energy) grid points that are dropped when calculating integral approximation with lower degree of precision. 



1. Model collision operator 

GS2 uses a model Fokker-Planck collision operator that includes the effects of pitch-angle scattering a nd en ergy 
diffusion while satsifying Boltzmann's H-Thcorem and conserving particle number, momentum, and energjEEJl 



where 



is the Lorentz collision operator, 



C[h] = C[h]+V[h] + M[h], (15) 



v d ( d /1 2X dh t l d 2 h 



Ax 2 dx \ dx Fq 

is the energy diffusion operator, and contains momentum- and energy-conserving corrections. The velocity- 

dependent collision frequencies v s and vn are given by 



ErfN - I (18) 
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and 



Erf[x] 



(19) 



with v a p the frequency of collisions of particles of species a with particles of species (3. A detailed description of the 
collision operator is given in Ref. 11 Here we simply present the gyroaveraged collision operator in spectral form: 
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where A: is the perpendicular wavenumber, a = fcuj_/Oo, Ai/ = i/p — v s , and 



+ Ji{a)v ± '- 



2v, 



a/3 
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4xe 



(21) 



Details on numerical implementation of the collision operator ( 20 1 can be found in Ref. [T2"l 



2. Numerical dissipation 



Numerical dissipation enters in GS2 through two mechanisms. The first is the optional decentering of spatial and 
temporal finite differences, as described in Ref. ITU] The lowest order contribution to dissipation due to decentering 
in time and space is 



d 2 h 

dm 



A0 15- 



1 



d 2 (x) 



R 



dtdO 



A6 [5- 



qFh 
T 



(22) 



where A9 is the grid spacing along the field line, At is the time step size, 5 is a parameter that allows for spatial 
upwinding (when 5 ^ 1/2), and (3 is a parameter that allows for the variation of the time discretization between fully 
explicit (f3 = 0) and fully implicit (/3 = 1)P. 

In order to see how this term leads to dissipation, we consider the simplified system governed by the equation 



Q. 



(23) 



dh dh 

Finite differencing this equation using the scheme given in Ref. 1101 we find that numerically we are solving the equation 



dh 
dt 



dh 
'd9 



d 2 h 
dtdO 



A9 5 



Assuming h = h(t)e ■ , we obtain the solution 



h(t) ~ exp 



kvt 



vAt[/3- 



1 



(24) 



(25) 



- k {AO (6 - 1/2) + vAt (P - 1/2)) _ 

which is damped unless (3 = 5= 1/2, as show in Fig. [5] While decentering of finite differences can sometimes improve 
numerical stability, care must be taken to ensure such artificial dissipation does not lead to unphysical behavior. This 
is typically done by monitoring the ratio of artificial to physical dissipation, which, ideally, should be small. 

The second source of numerical dissipation arises in systems with sheared magentic fields due to the necessity of 
a 'twist-and-shift' parallel boundary conditional This non-periodic boundary condition couples modes at opposite 
ends of the simulation domain along the field line. Since only a finite number of modes can be kept in a simulation, 
some modes will eventually couple to modes that are not present, and this information is lost. The information that 
is lost is replaced by a smoothed distribution function, which should be associated with an increase in the entropy of 
the system. This entropy generation should be diagnosed in order to verify that it is small compared to the entropy 
generated by collisions. 
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IV. VELOCITY SPACE RESOLUTION DIAGNOSTICS 



There are numerous ways in which one could try to determine whether or not a particular simulation is well-resolved 
in velocity space. Ideally, one would perform a grid convergence study for each simulation; if quantities of interest are 
unchanged by doubling the number of grid points, one can feel relatively confident in the simulation results. However, 
this process is computationally expensive, as it involves running a simulation multiple times with an excessive number 
of grid points. Consequently, it is not desirable to perform a grid convergence study for every simulation. In practice, 
one tests convergence for a problem thought to be resolution intensive and posits that other simulations, which likely 
require fewer grid points, are therefore resolved. Unfortunately, one seldom knows in advance how fine the structure 
in velocity space will become, so one can't be fully confident that every simulation is resolved. 

An alternative approach that has recently gained popularity in the computational plasma physics community in- 
volves monitoring entropy balance in the systerrpl2l. The entropy balance relation arises from multiplying the gyroki- 
netic equation (JTJ by HTq/Fq and integrating over all phase space. Since the gyrokinetic equation itself is automatically 
satisfied by a gyrokinetic solver, the only possible sources of inbalance in this relation come from numerical dissipa- 
tion and errors in the numerical approximations to phase space integrals. If the change in entropy due to numerical 
dissipation is also diagnosed and included in the entropy balance, as is often the case, then we are left with errors due 
only to phase space integration. Since the errors in these particular integrals are not directly related to errors in the 
calculation of the distribution function at the newest timestep, they do not necessarily correlate with the simulation 
resolution. In particular, one could easily define a poorly-resolved system for which this diagnostic predicts perfect 
entropy balance. One such example is the linear, collisionless ion acoustic wave in a slab (treated in detail in Appendix 
A). For this case, we numerically find entropy balance despite the fact that the numerical damping rate goes bad due 
to poor resolution in velocity space. 
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Of course, one could simply produce plots or movies of the distribution function in velocity space over the course of 
the simulation to see if structure develops at the gridscale. This is undoubtedly useful and possibly sufficient in some 
cases. However, what exactly one sees depends on how the data is visualized; for data on irregularly spaced grids, the 
interpolation scheme used to generate the images often introduces erroneous or misleading structure. Furthermore, 
for simulations involving non-trivial spatial structure, one would have to examine movies of the distribution function 
at each point in physical space. This is a memory- and time-intensive approach that is rarely feasible. 

We would like to have computationally cheap diagnostics that provide real-time information on velocity space 
resolution that is easy to analyze and interpret. In the following subsections, we present two such diagnostics developed 
for implementation in GS2 that could easily be adapted for use in other continuum kinetic simulations. 



A. Integral error estimates 



Upon consideration of the collisionless gyrokinetic-Maxwell's system of equations, one finds that the only nontrivial 
operation in velocity space is integration, which enters in the calculation of the electromagnetic fields. Consequently, 
resolution in velocity space is limited only by the accuracy with which the velocity space integrals are calculated. By 
calculating the error in our numerical integration, we are thus able to monitor velocity space resolution. 

In particular, when we discretize the gyrokinetic equation, we obtain an equation of the form 

9j+i = G [9j , $j , $j+u Xj , Xj+i] > (26) 

where g = (/i) is the perturbed, guiding center distribution function evolved by GS2, <!> is the electrostatic potential, 
X is the generalized electromagnetic potential defined in Eqn. ([5]), G is a function that depends on the details of 
the numerical scheme, and the subscript denotes the timestep. We assume that the time-converged solution for g is 
independent of the initial condition. Since using the calculated <?j and <frj is equivalent to specifying a new initial 
condition, we find that the time-converged solution is independent of errors in g and $ at earlier timesteps. This is 
convenient because it means we can monitor resolution merely by calculating the error made in the latest timesteps 
of a time-converged simulation. 

Ideally, we would accomplish this by calculating estimates for the error in and Xj+i an d plugging these into 



Eqn. (26) to obtain an error estimate for gj+i- This might be feasible for linear systems, but the presence of nonlinear 
terms makes this approach computationally prohibitive. Consequently, we must define an alternative quantity whose 
error estimate is cheaper to compute, but that can still be used as a means of monitoring velocity space resolution. 
There are numerous possible candidates; we choose to compute two quantities, i>$ and va, related to Vj^4> and Vj_^4y: 

where k x and ky are the wavenumbers corresponding to the coordinates x = (tp — ipo) qa/B^r^ and y = 
— (a — a ) ro/qo 1 ®- Here, is the poloidal flux, a is the field line label, B is the background magnetic field at 
the magnetic axis, tq is the distance from the magnetic axis to the center of the simulation domain, and q$ is the 



safety factor on the field line of interest, labeled by (ipo, ot^). The quantities in Eqn. (27l were chosen because, with the 
exceptions of the parallel convection term and one source term, 4> and A\\ always enter the gyrokinetic equation for g 
multiplied by either k x or k y . Therefore, it is reasonable that this k- weighted quantity is most likely to be responsible 



for errors in gj+i- Although not considered here, the expression (27 1 could potentially be improved by including 
in the max operator. This would take into account the effect of the parallel convection term. However, recent 
theoretical and numerical work suggests that velocity space structure may be generated primarily by nonlinear 
perpendicular phase mixing (instead of linear, parallel phase mixing). 

Having chosen appropriate indicators of velocity space resolution, we must devise a method for estimating the error 
in these quantities. This error depends on the particular numerical integration scheme used. For the energy and 
untrapped A integrals, which use Gaussian quadrature, the error, cq, is given by 

e G = 7 m / (2m) (C), (28) 

where / is the integrand, m is the number of grid points, and £ is some unknown point in the interval of integration. 
The quantity 7 m is 
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for the untrapped A and finite domain energy integrals that use Gauss-Legendre quadrature and 

(ml) 2 



(2m)! 



(30) 



for the semi-infinite domain energy integral that uses Gauss-Laguerre quadrature. The error, cl, for the trapped 
A integrals, which use a newly upgraded integration scheme based on Lagrange interpolating polynomials (see e.g. 
Ref. [16]), is given by 



where 



1 f f(m) {07T{z)dZj 

ml J 

m 

K(z) = Y[(Z Zi) , 



(31) 



(32) 



with Zi the i th grid point. It should be noted that £ in Eqn ( |3l| is an unknown function of z whose domain is some 
subset of the interval of integration. 

From Eqns (28 1 and (31 1, we see that Gaussian quadrature gives exact results for polynomials of degree less than 



2m, while the Lagrangian method gives exact results only for polynomials of degree less than m. We say that the 
two schemes have degrees of precision 2m — 1 and m — 1, respectively. This difference arises because the grid points 
in the Lagrangian method are fixed by boundary conditions, whereas the grid points in Gaussian quadrature are free 
parameters optimally chosen to improve the scheme's degree of precision. 

Unfortunately, the formal error expressions (28) and (31 1 are not very useful in practice: they require information 



about high-order derivatives of the distribution function, which is unavailable. As an alternative estimate for the 
error, we choose to compare multiple integral approximations computed with different degrees of precision, a common 
technique in numerical analysis^. 



1. General description of the scheme 

Given the value of a function f(z) at N fixed points on the interval [a, b], we would like to find two different 
approximations to the integral J f(z)dx. In our earlier discussion, we stated that an approximation with degree of 
precision N — 1 can be found using a technique based on Lagrange interpolation; we call this approximation Ah- If we 
instead choose to use only M of the given functional values (M < N) , we can use the same technique to find another 
integral approximation, Ai, with degree of precision M — 1. An estimate for the absolute error e a in the less accurate 
of these two approximations is obtained by taking the difference between the two: 

€ a = \A h -A l \. (33) 

Making the reasonable assumption that the approximation with higher degree of precision is more accurate, e a 
represents the error in Ai. However, it can also be used as a more conservative error estimate for Ah- 

If the N points are chosen according to Gaussian quadrature rules, then one can find an integral approximation 
with degree of precision 27V — 1. As before, a second approximation can be obtained by using only M of the N grid 
points. However, due to the uniqueness of the grid points used for Gaussian quadrature, the M-point grid no longer 
satisfies Gaussian quadrature rules. As a result, this second approximation once again has degree of precision M — 1. 
Since the degrees of precision of the two approximations differ by greater than a factor of two, the resulting error 
estimate is likely to be very conservative when applied to Ah- The factor of approximately two difference in degree 
of precision makes this error estimate similar to that obtained by comparing results from runs with N and N/2 grid 
points, respectively (for which the degrees of precision would be 2N — 1 and TV — 1). 

The conservative nature of the error estimate for Ah depends upon our assumption that a higher degree of precision 
results in a more accurate integral approximation. For Gaussian quadrature, it can be shown that the error in the 
integral approximation can be made arbitrarily small by choosing the degree of precision large enough^. The same 
result does not necessarily hold for the Lagrangian method with arbitrary grid spacing because the weights in this 
case are not all guaranteed to be positive. However, the error cm in an M-point integral approximation satisfies 

M 

e« <2 e £|t«,M 

i=l 

< 2eM max \wA 
= 2cMk{M), 



(34) 

(35) 
(36) 



12 



where e can be chosen arbitrarily small for large enough M, and w\ is the weight corresponding to the i grid 
point out of M. From this result, we see that as long as k is bounded when M — ► oo, then em -> as M -> oo. 
This cannot be verified in advance, but one can gain confidence by checking a posteriori. In practice, we calculate k 
for the chosen M and subdivide the integration domain into subintervals with fewer points if k is larger than some 
reasonable value. 



2. Implementation in GS2 

In GS2, we must compute two-dimensional integrals over energy and A. As stated in Sec. Ill, each of these integrals 
is effectively separated into two by splitting the A integration into trapped and untrapped regions. Since the number 
of grid points in energy and both A regions can be varied independently of each other, we wish to monitor resolution 
in each of these three variables individually. This entails computing three separate integral error estimates: one for 
energy integrals, one for untrapped A integrals, and one for trapped A integrals. 

These integral error estimates are calculated using the technique described in the previous subsection. For energy 
and untrapped A integrals, Gaussian quadrature is used to obtain the two-dimensional integral approximation Ah- 
This approximation has degree of precision 2Ne — 1 for the energy integration and 2N U — 1 for the untrapped A 
integration, where Ne and N u are the number of energy and untrapped A grid points, respectively. To obtain the 
second approximation, Ai, we fix the grid and weights for one variable and drop one grid point for the other variable, 
recomputing the weights. As an example, we choose to drop an untrapped A grid point. The degree of precision for A[ 
is then 2Ne ~ 1 for the energy integration and N u — 2 for the untrapped A integration. Since there is nothing special 
about the particular grid point we drop, we repeat the process a total of N u times, each time dropping a different 
point and computing a different set of weights. The final error estimate is an average of these error estimates. 

For the trapped A integrals, Lagrangian quadrature is used to obtain Ah, which has degree of precision JV* — 1. 
We obtain the approximation Ai by dropping two points symmetrically about v\\ = 0, as shown in Fig. |4j We drop 
an additional point here because it provides a slightly more conservative error estimate and because maintaining 
the symmetry of the grid points provides better stability for the weights associated with the Lagrange interpolation 
scheme. As before, we repeat this process for each possible grid point pair and take the average of the individual error 
estimates to get the final error estimate. 

All modified grids and weights necessary for the integral error estimates are computed once at initialization and 
need not be computed again. The additional integrations necessary to obtain our error estimates are computationally 
cheap when compared to the expense of solving for the distribution function and fields at each time step. Furthermore, 
we do not need an error estimate at each time step, so the diagnostic can be used sparingly. Consequently, our error 
estimate comes at essentially no extra cost. 



B. Spectral method 



An alternative method for testing v-space resolution is to expand the velocity space distribution function in an 
appropriate basis set and monitor the amplitude of the basis function coefficients. Whenever the highest mode 
number coefficients that can be accurately calculated in the simulation acquire appreciable amplitudes, we can no 
longer feel confident that the simulation is resolved. Since we choose our grid points according to Gauss-Legendre 
quadrature, it is convenient (and most accurate) to choose the Legendre polynomials as our basis functions. The 
coefficient of the m th Legendre polynomial in the expansion of h is given by 



2m + 1 
2 

2m + 1 
2 



l 

h{s)P m (s)ds (37) 



i 

n.—l 



w.Ms t )P m (sA, (38) 



i=l 



where P m is the m Legendre polynomial, and {wi} are the weights associated with Gauss-Legendre quadrature. The 



integral approximation in Eqn. (38 1 has degree of precision 2N — 1. Assuming h has a degree of at least m (otherwise 
c m = 0), our approximation for c m is only exact for m < N. 

There are various ways in which one could use these {c m } to estimate the error in velocity space resolution. We 
assume locality of interaction between the various modes so that we only have to monitor the amplitudes of the few 
highest modes. At each (0, k x , fc y )-point, we find the maximum amplitude of the three highest mode number spectral 
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FIG. 6: Comparison of actual and (unsealed) estimated error in wave frequency due to insufficient resolution in energy (left) 
and untrapped A (right). The actual wave frequency, w, is determined from a higher resolution run with 64 grid points in 

energy and both trapped and untrapped A. The actual relative error, e, is then defined to be e = <\J where ui n is the 

approximation to uo obtained from a run with n grid points. 



coefficients, Ch.max, and the maximum amplitude of all the spectral coefficients, c max . We then use the following 
normalized sum as a relative estimate for the error: 

^ ^ Ch,niax(@ ) kx> ^y) / ^ ^ Cmax(@ : kx: ky) • (39) 

When the normalized amplitude e c grows too large, we can no longer be confident that the simulation is resolved. 
Of course, how large e c can get before resolution suffers varies from problem to problem. As before with the integral 
method, we determine a scaled estimate of the error based on empirical evidence from a wide range of simulation 
data. 



1. Application of error diagnostics 



We have applied both the integral and spectral error diagnostics to a diverse set of simulations, including: linearly 
growing modes such as the electron drift wave and the ITG mode; linearly damped modes such as the ion acoustic 
wave and kinetic Alfven wave; neoclassical transport; and nonlinear dynamics of slab ETG and toroidal ITG modes. 
From these simulations, we have determined empirical scaling factors for our conservative error estimates. Here, we 
present typical results from a cross-section of the above simulations. 

Fig. [6] compares the unsealed error estimates in energy and A with the actual errors in growth rate as we vary 
the number of grid points in a linear simulation of the collisionless toroidal ITG mode (using Cyclone base case 
parameters^). The simulation remains well-resolved down to very few grid points, and the error estimates agree well 
with the actual error. The error due to resolution in untrapped A is still small for as little as four grid points due to 
our choice of velocity variables, as illustrated by the snapshot of the distribution function shown in Fig. [7] 

Figs.[H]and[9]show the damping of Ay and the corresponding scaled error estimates for the simulation of a collisionless 
kinetic Alfven wave with 16 energy grid points and 32 pitch angles for each sign of the parallel velocity. The collisionless 
damping rate in Fig. [8] agrees with theory until sub-gridscale structure develops in velocity space, at which point 
damping ceases. The onset of sub-gridscale structure corresponds to the peak in scaled error in Fig. [9] The addition 
of a small collisionality prevents sub-gridscale structure, as shown in Fig. [8] where the damping rate of Ay agrees well 
with theory indefinitely. This is accurately predicted by the error estimates of Fig. |10| which never reach appreciable 
magnitude. 
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FIG. 7: Non-Boltzmann part of the perturbed distribution function, normalized by Jo, for the linear, toroidal ITG mode with 
Cyclone base-case parameters. The use of a polar grid in velocity space, as well as the fine mesh near the trapped-passing 
boundary, minimizes the number of grid points necessary for resolution. 

V. ADAPTIVE COLLISION FREQUENCY 

As stated earlier, we would like to know what combination of dissipation and grid spacing is necessary for a resolved 
simulation. One way to approach this problem is to fix the dissipation and vary the number of grid points to find 
how many are required to get an accurate result. This is the general idea behind the error estimation diagnostics 
described in the previous section. However, if we wanted to use this approach to ensure that the simulation remained 
resolved, we would have to implement an adaptive grid, which is difficult to do on massive, multi-processor machines. 

Instead, we choose an alternative approach: we fix the number of grid points and vary the dissipation until we 
have a well-resolved result. In particular, we have implemented an adaptive collision frequency in GS2 that allows for 
the independent variation of the collisionality associated with pitch-angle scattering and energy diffusion. Given an 
acceptable error tolerance for velocity space calculations, a scaled version of the integral error estimate described in 
the previous section is used to determine whether or not the simulation is well-resolved. The collision frequency is 
then adjusted using a feedback process until the scaled estimate of the error converges to within some pre-specified 
window of the desired error tolerance. In this way, the approximate minimum possible dissipation is used to achieve 
an acceptable degree of resolution in velocity space. 

Of course, the amount of dissipation necessary to resolve a simulation at a fixed number of grid points may be quite 
large if a coarse grid is used. Consequently, the collisionless dynamics may be modified. As a result, it is necessary 
to compare the converged collision frequency with dynamic frequencies of interest in the problem. 

As an example we consider a nonlinear simulation of electron temperature gradient (ETG) turbulence in slab 
geometry (i.e. straight background magnetic field). In the nonlinear phase, small scales are expected to develop in 
velocity space, potentially challenging numerical resolution. In Fig.[TT] we see that this is indeed the case. Our velocity 
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Time [(k„v lhii )-'] Time [( k iVi)~'] 

FIG. 8: Barnes damping of the kinetic Alfven wave for simulations with 16 energy grid points and 32 pitch angles for each 
sign of In the absence of collisions (left), sub-grid scale structures develop in velocity space, and the damping is artificially 
terminated. A small collisionality (y <C 7) prevents the development of sub-grid scale structures in velocity space, and the 
damping rate remains correct indefinitely (right). 
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FIG. 9: Time evolution of the integral (left) and spectral (right) error estimates for the collisionless kinetic Alfven wave. Vertical 
line represents time at which damping rate artificially terminates due to poor resolution. 



space resolution diagnostics indicate that the errors in velocity space begin to increase sharply during the transition 
from linear instability to turbulence. However, our use of an adaptive collision frequency prevents the estimated error 
from exceeding the user-defined relative error tolerance (in this case, 0.01). We see that the error remains on the 
threshold of the error tolerance, while the collision frequency for energy diffusion increases to a steady-state value of 
v w 0.027 k\\Vth,e! which is well below the dynamic frequency in the system. Consequently, the collisionless dynamics 
are unaltered. 
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FIG. 10: Time evolution of integral and spectral error estimates for the weakly collisional kinetic Alfven wave damping. The 
estimates correctly indicate that the simulation remains well resolved indefinitely. 
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FIG. 11: (Left): Normalized electron heat flux vs. time for a nonlinear simulation of ETG turbulence. Scaled estimates of the 
error in energy and A resolution increase during nonlinear saturation, but are kept within the specified relative error tolerance 
of 0.01 with the use of an adaptive collision frequency. (Right): Collision frequency (normalized by k\\Vth,e) vs. time. 



VI. SUMMARY 



In this paper, we discussed the development of small-scale structure in velocity space, presented a set of velocity 
space resolution diagnostics for use in gyrokinetic simulations, and introduced an adaptive collisionality that allows 
us to resolve simulations with an approximate minimal necessary dissipation for a fixed number of grid points in 
velocity space. In Sec. |ll]we demonstrated the tendency of collisionless plasmas to develop increasingly fine scales in 
the distribution of particle velocities and discussed the phase mixing processes that lead to such behavior. 

In Sec. |TTT] we described the treatment of velocity space in the gyrokinetic code GS2. We gave details on the 
choice of velocity space variables (energy and pitch-angle) and discretization scheme, which is chosen to minimize the 
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error of the numerical integrals necessary to obtain the electromagnetic fields. This included presentation of a newly 
implemented energy grid, which provides spectrally accurate integrals over particle energies. Additionally, we gave a 
brief discussion of both the physical and numerical dissipation mechanisms available for use in GS2. 

We discussed common approaches to monitoring velocity space resolution in Sec. |IV| and the difficulties associated 
with each. We then proposed two new measures of velocity space resolution and detailed implementation in GS2. 
One of the proposed resolution diagnostics involves obtaining estimates for the error in field integrals by comparing 
numerical integrals obtained using integration schemes with differing degrees of precision. The other resolution 
diagnostic involves decomposing the perturbed distribution function into spectral components in velocity space and 
monitoring the amplitude of the spectral coefficients. Both diagnostics should be quite conservative. 

We then applied our resolution diagnostics to a number of example problems, including Landau damping of the ion 
acoustic wave, Barnes damping of the kinetic Alfven wave, and linear instability of the toroidal ITG mode. We found 
that both diagnostics do well in qualitatively estimating errors due to limited velocity space resolution. Due to their 
conservative nature, an empirical scaling factor was necessary to obtain correct quantitative predictions. 

In Sec. [V] we coupled the error estimates from our resolution diagnostics with a model physical collision operator to 
develop an adaptive collision frequency. This adaptive collision frequency allowed us to resolve velocity space while 
using an approximate minimal necessary amount of dissipation. When using the adaptive collision frequency, one 
must monitor the ratio of the collision frequency to the dynamic frequency to ensure that one is still within the weakly 
collisional regime. 

In conclusion, we found that dissipation was not necessary to resolve linear instabilities, but it was necessary to 
resolve nonlinear dynamics and linearly damped waves. For the nonlinear cases considered here (slab ETG and 
toroidal ITG), the required collisionality for resolution obtained with the adaptive collision frequency was found to 
be no larger than the physical collisionality used in modern fusion experiments. 



APPENDIX A: LANDAU-DAMPED ION ACOUSTIC WAVE 



We consider the collisionless ion acoustic wave in slab geometry with adiabatic electrons. The gyrokinetic equation 
for this system has the particularly simple form 

Changing variables from h to g = (/i) and assuming solutions of the form 

g=~g(v)e i ( k "*-" t ), (A2) 

we obtain 

(u-kv)g = kv e -^F M , (A3) 
where we are using v = v\\ and k — k\\ for convenience. Neglecting FLR effects and assuming quasineutrality gives 

(u -kv)g = kvr— [ d 3 v'g(v'). (A4) 
n J 

Defining 



g{v) =2n vj_dvj_g(v) (A5) 
Jo 

and integrating over the perpendicular velocities in the gyrokinetic equation yields 

(w-kv)g{v) = kF(v), (A6) 

where 

F(v) = vt "l 1 e~^, (A7) 



m = / dv'g(v'). (A8) 
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Following the analysis of Refs. [55] and HS1 we see that this equation has solutions of the form 

1 



g(v) = F(v) 



u — v 



+ A(k,u)S(u - v) 



with u = j, provided that A is chosen to satisfy the condition 

m = [ dv'g(v')=V [ dv'^^- + A(k,u)F(u). 
J J u — v 

A general solution is given in the form 

/OO POO 
/ C(k,u)g ktU (v)e ik ^- ut Uk du, 
-co J — OO 

where C(k,u) is determined by the initial condtion 

J(z,v,Q) = J J C(k,u)g k Jvy kz dk du. 
Taking the inverse Fourier transform of the above expression gives 

F(k,v) = J C(k,u)g ktU (v)du, 

where 

F(k,v) = ~ f~f(z,v,Q)e- lkz dz. 



Plugging the expression (A9) for g into the initial condition ( A13 1 yields 



F(v) 



= V l ^hA du + k(k,v)C{v). 



(A9) 



(A10) 



(All) 



(A12) 



(A13) 



(A14) 



(A15) 



We now have two equations, (AlOl and (A15l, for two unknowns (A and C). In order to solve this linear system, 
it is convenient to define some new notation. Any square integrable function H can be written 



/oo 
K(p)e lpq dp. 
-OO 



We define the positive and negative frequency parts of H as 

f ±oo 



P it oo 

H±(q) = ± / K( P yPidp, 
Jo 



(A16) 



(A17) 



so that H = H + + H- . Further we define the function = H + — H- . It can be shown that has the alternate 
form 



m J~oc V -V 



With these definitions in hand, we rewrite eqns (AlOl and (A15) as 



m = -niF*(u) + AF(u), 
T{k,v) 



F(v) 



= (A + TTi)C+(v) + (A-TTi)C-(v). 



Eliminating A gives an expression involving C+ and C_ : 

F(k,u) = (ni +2mF + (u))C+{k,u) + (m - 2iriF^(u)) C_(fc, u). 



(A18) 

(A19) 
(A20) 

(A21) 
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The transform T can also be broken down into negative and positive frequency parts to give two separate equations. 



F±{k,u) = {n x ±2mF±{u))C±{u) 
These can then be used to construct C(k, u): 

u) T-{k,u) 



C(k,u) 



ni + 2itiF + (u) rii — 2niF-(u) 



Substituting the expressions {A9\ and {A23\ for g and C into the equation (All) for / gives 



f(z,v,t) 



ni + 2niF + (u) m — 2mF-(u) 
1 



F(v) 



V \-A(k,u)5(u-v) 



ik(z-ut) dk du _ 



We can use the identity 



T±(k,u) 



1 

2n 



-ikz' jiJ 



dz' / S±(u-i/)f(z',v',0)dv' 



to rewrite eqn ( A24 1 in the more convenient form 

6-\-(u — v') 5-(u — v') 



f(z,v,t) 



n\ + 2'KiF + (u) m — 2iriF-(u) 



f(z',v',0) 
2tt 



F(v) 



1 



+ A(fe, u)5(u - v) 



e ik ( z - z '- ut )dz'dv'dk du. 



Now we pick an initial condition of the form 



which gives 



f(z,v,0) = f(v,0)e ik ° z , 



7(z,v,t) = e»*0-*> [m+viFM) I /+( "' 0) , , + J 
Jy ' ' ' y \n 1 + 2niF + (v) n x - 2mF-{v) J 

m ( + /-( U ,o) ) jM*-^ du . 

u — v \ m + 2iriF + (u) n\ — 2niF-(u) J 



(A22) 



(A23) 



(A24) 



(A25) 



(A26) 



(A27) 

(A28) 
(A29) 
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